home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 15 / CU Amiga Magazine's Super CD-ROM 15 (1997)(EMAP Images)(GB)[!][issue 1997-10].iso / CUCD / Graphics / Ghostscript / source / gsmisc.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-07-01  |  23.8 KB  |  818 lines
  1. /* Copyright (C) 1989, 1996, 1997 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of Aladdin Ghostscript.
  4.   
  5.   Aladdin Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author
  6.   or distributor accepts any responsibility for the consequences of using it,
  7.   or for whether it serves any particular purpose or works at all, unless he
  8.   or she says so in writing.  Refer to the Aladdin Ghostscript Free Public
  9.   License (the "License") for full details.
  10.   
  11.   Every copy of Aladdin Ghostscript must include a copy of the License,
  12.   normally in a plain ASCII text file named PUBLIC.  The License grants you
  13.   the right to copy, modify and redistribute Aladdin Ghostscript, but only
  14.   under certain conditions described in the License.  Among other things, the
  15.   License requires that the copyright notice and this notice be preserved on
  16.   all copies.
  17. */
  18.  
  19. /* gsmisc.c */
  20. /* Miscellaneous utilities for Ghostscript library */
  21. #include "malloc_.h"
  22. #include "math_.h"
  23. #include "memory_.h"
  24. #include "gx.h"
  25. #include "gpcheck.h"        /* for gs_return_check_interrupt */
  26. #include "gserrors.h"
  27. #include "gconfigv.h"        /* for USE_ASM */
  28. #include "gxfarith.h"
  29. #include "gxfixed.h"
  30.  
  31. /* Define private replacements for stdin, stdout, and stderr. */
  32. FILE *gs_stdin, *gs_stdout, *gs_stderr;
  33. /* Ghostscript writes debugging output to gs_debug_out. */
  34. /* We define gs_debug and gs_debug_out even if DEBUG isn't defined, */
  35. /* so that we can compile individual modules with DEBUG set. */
  36. char gs_debug[128];
  37. FILE *gs_debug_out;
  38.  
  39. /* Define eprintf_program_name and lprintf_file_and_line as procedures */
  40. /* so one can set breakpoints on them. */
  41. void
  42. eprintf_program_name(FILE *f, const char *program_name)
  43. {    fprintf(f, "%s: ", program_name);
  44. }
  45. void
  46. lprintf_file_and_line(FILE *f, const char *file, int line)
  47. {    fprintf(f, "%s(%d): ", file, line);
  48. }
  49.  
  50. /* Log an error return.  We always include this, in case other */
  51. /* modules were compiled with DEBUG set. */
  52. #undef gs_log_error        /* in case DEBUG isn't set */
  53. int
  54. gs_log_error(int err, const char _ds *file, int line)
  55. {    if ( gs_log_errors )
  56.       { if ( file == NULL )
  57.           dprintf1("Returning error %d.\n", err);
  58.         else
  59.           dprintf3("%s(%d): Returning error %d.\n",
  60.                (const char *)file, line, err);
  61.       }
  62.     return err;
  63. }
  64.  
  65. /* Check for interrupts before a return. */
  66. int
  67. gs_return_check_interrupt(int code)
  68. {    if ( code < 0 )
  69.       return code;
  70.     { int icode = gp_check_interrupts();
  71.       return (icode == 0 ? 0 :
  72.           gs_note_error((icode > 0 ? gs_error_interrupt : icode)));
  73.     }
  74. }
  75.  
  76. /* ------ Substitutes for missing C library functions ------ */
  77.  
  78. #ifdef memory__need_memmove        /* see memory_.h */
  79. /* Copy bytes like memcpy, guaranteed to handle overlap correctly. */
  80. /* ANSI C defines the returned value as being the src argument, */
  81. /* but with the const restriction removed! */
  82. void *
  83. gs_memmove(void *dest, const void *src, size_t len)
  84. {    if ( !len )
  85.         return (void *)src;
  86. #define bdest ((byte *)dest)
  87. #define bsrc ((const byte *)src)
  88.     /* We use len-1 for comparisons because adding len */
  89.     /* might produce an offset overflow on segmented systems. */
  90.     if ( ptr_le(bdest, bsrc) )
  91.       {    register byte *end = bdest + (len - 1);
  92.         if ( ptr_le(bsrc, end) )
  93.           {    /* Source overlaps destination from above. */
  94.             register const byte *from = bsrc;
  95.             register byte *to = bdest;
  96.             for ( ; ; )
  97.               {    *to = *from;
  98.                 if ( to >= end )    /* faster than = */
  99.                   return (void *)src;
  100.                 to++; from++;
  101.               }
  102.           }
  103.       }
  104.     else
  105.       {    register const byte *from = bsrc + (len - 1);
  106.         if ( ptr_le(bdest, from) )
  107.           {    /* Source overlaps destination from below. */
  108.             register const byte *end = bsrc;
  109.             register byte *to = bdest + (len - 1);
  110.             for ( ; ; )
  111.               {    *to = *from;
  112.                 if ( from <= end )    /* faster than = */
  113.                   return (void *)src;
  114.                 to--; from--;
  115.               }
  116.           }
  117.       }
  118. #undef bdest
  119. #undef bsrc
  120.     /* No overlap, it's safe to use memcpy. */
  121.     memcpy(dest, src, len);
  122.     return (void *)src;
  123. }
  124. #endif
  125.  
  126. #ifdef memory__need_memchr        /* see memory_.h */
  127. /* ch should obviously be char rather than int, */
  128. /* but the ANSI standard declaration uses int. */
  129. const char *
  130. gs_memchr(const char *ptr, int ch, size_t len)
  131. {    if ( len > 0 )
  132.     {    register const char *p = ptr;
  133.         register uint count = len;
  134.         do { if ( *p == (char)ch ) return p; p++; } while ( --count );
  135.     }
  136.     return 0;
  137. }
  138. #endif
  139.  
  140. #ifdef memory__need_memset        /* see memory_.h */
  141. /* ch should obviously be char rather than int, */
  142. /* but the ANSI standard declaration uses int. */
  143. void *
  144. gs_memset(void *dest, register int ch, size_t len)
  145. {    if ( ch == 0 )
  146.         bzero(dest, len);
  147.     else if ( len > 0 )
  148.     {    register char *p = dest;
  149.         register uint count = len;
  150.         do { *p++ = (char)ch; } while ( --count );
  151.     }
  152.     return dest;
  153. }
  154. #endif
  155.  
  156. #ifdef malloc__need_realloc        /* see malloc_.h */
  157. /* Some systems have non-working implementations of realloc. */
  158. void *
  159. gs_realloc(void *old_ptr, size_t old_size, size_t new_size)
  160. {    void *new_ptr;
  161.     if ( new_size )
  162.       { new_ptr = malloc(new_size);
  163.         if ( new_ptr == NULL )
  164.           return NULL;
  165.       }
  166.     else
  167.       new_ptr = NULL;
  168.     /* We have to pass in the old size, since we have no way to */
  169.     /* determine it otherwise. */
  170.     if ( old_ptr != NULL )
  171.       { if ( new_ptr != NULL )
  172.           memcpy(new_ptr, old_ptr, min(old_size, new_size));
  173.         free(old_ptr);
  174.       }
  175.     return new_ptr;
  176. }
  177. #endif
  178.  
  179. /* ------ Debugging support ------ */
  180.  
  181. /* Dump a region of memory. */
  182. void
  183. debug_dump_bytes(const byte *from, const byte *to, const char *msg)
  184. {    const byte *p = from;
  185.     if ( from < to && msg )
  186.         dprintf1("%s:\n", msg);
  187.     while ( p != to )
  188.        {    const byte *q = min(p + 16, to);
  189.         dprintf1("0x%lx:", (ulong)p);
  190.         while ( p != q ) dprintf1(" %02x", *p++);
  191.         dputc('\n');
  192.        }
  193. }
  194.  
  195. /* Dump a bitmap. */
  196. void
  197. debug_dump_bitmap(const byte *bits, uint raster, uint height, const char *msg)
  198. {    uint y;
  199.     const byte *data = bits;
  200.     for ( y = 0; y < height; ++y, data += raster )
  201.         debug_dump_bytes(data, data + raster, (y == 0 ? msg : NULL));
  202. }
  203.  
  204. /* Print a string. */
  205. void
  206. debug_print_string(const byte *chrs, uint len)
  207. {    uint i;
  208.     for ( i = 0; i < len; i++ )
  209.       dputc(chrs[i]);
  210.     fflush(dstderr);
  211. }
  212.  
  213. /* ------ Arithmetic ------ */
  214.  
  215. /* Compute M modulo N.  Requires N > 0; guarantees 0 <= imod(M,N) < N, */
  216. /* regardless of the whims of the % operator for negative operands. */
  217. int
  218. imod(int m, int n)
  219. {    if ( n <= 0 )
  220.       return 0;        /* sanity check */
  221.     if ( m >= 0 )
  222.       return m % n;
  223.     { int r = -m % n;
  224.       return (r == 0 ? 0 : n - r);
  225.     }
  226. }
  227.  
  228. /* Compute the GCD of two integers. */
  229. int
  230. igcd(int x, int y)
  231. {    int c = x, d = y;
  232.     if ( c < 0 ) c = -c;
  233.     if ( d < 0 ) d = -d;
  234.     while ( c != 0 && d != 0 )
  235.       if ( c > d ) c %= d;
  236.       else d %= c;
  237.     return d + c;        /* at most one is non-zero */
  238. }
  239.  
  240. #if defined(set_fmul2fixed_vars) && !USE_ASM
  241.  
  242. /*
  243.  * Floating multiply with fixed result, for avoiding floating point in
  244.  * common coordinate transformations.  Assumes IEEE representation,
  245.  * 16-bit short, 32-bit long.  Optimized for the case where the first
  246.  * operand has no more than 16 mantissa bits, e.g., where it is a user space
  247.  * coordinate (which are often integers).
  248.  *
  249.  * The assembly language version of this code is actually faster than
  250.  * the FPU, if the code is compiled with FPU_TYPE=0 (which requires taking
  251.  * a trap on every FPU operation).  If there is no FPU, the assembly
  252.  * language version of this code is over 10 times as fast as the emulated FPU.
  253.  */
  254. /* Some of the following code has been tweaked for the Borland 16-bit */
  255. /* compiler.  The tweaks do not change the algorithms. */
  256. #if arch_ints_are_short && !defined(FOR80386)
  257. #  define SHORT_ARITH
  258. #endif
  259. int
  260. set_fmul2fixed_(fixed *pr, long /*float*/ a, long /*float*/ b)
  261. {
  262. #ifdef SHORT_ARITH
  263. #  define long_rsh8_ushort(x)\
  264.     (((ushort)(x) >> 8) | ((ushort)((ulong)(x) >> 16) << 8))
  265. #  define utemp ushort
  266. #else
  267. #  define long_rsh8_ushort(x) ((ushort)((x) >> 8))
  268. #  define utemp ulong
  269. #endif
  270.     /* utemp may be either ushort or ulong.  This is OK because */
  271.     /* we only use ma and mb in multiplications involving */
  272.     /* a long or ulong operand. */
  273.     utemp ma = long_rsh8_ushort(a) | 0x8000;
  274.     utemp mb = long_rsh8_ushort(b) | 0x8000;
  275.     int e = 260 + _fixed_shift - ((
  276.         (((uint)((ulong)a >> 16)) & 0x7f80) +
  277.         (((uint)((ulong)b >> 16)) & 0x7f80)
  278.       ) >> 7);
  279.     ulong p1 = ma * (b & 0xff);
  280.     ulong p = (ulong)ma * mb;
  281. #define p_bits (size_of(p) * 8)
  282.  
  283.     if ( (byte)a )        /* >16 mantissa bits */
  284.     {    ulong p2 = (a & 0xff) * mb;
  285.         p += ((((uint)(byte)a * (uint)(byte)b) >> 8) + p1 + p2) >> 8;
  286.     }
  287.     else
  288.         p += p1 >> 8;
  289.     if ( (uint)e < p_bits )        /* e = -1 is possible */
  290.         p >>= e;
  291.     else if ( e >= p_bits )        /* also detects a=0 or b=0 */
  292.       {    *pr = fixed_0;
  293.         return 0;
  294.       }
  295.     else if ( e >= -(p_bits - 1) || p >= 1L << (p_bits - 1 + e) )
  296.         return_error(gs_error_limitcheck);
  297.     else
  298.         p <<= -e;
  299.     *pr = ((a ^ b) < 0 ? -p : p);
  300.     return 0;
  301. }
  302. int
  303. set_dfmul2fixed_(fixed *pr, ulong /*double lo*/ xalo, long /*float*/ b, long /*double hi*/ xahi)
  304. {
  305. #ifdef SHORT_ARITH
  306. #  define long_lsh3(x) ((((x) << 1) << 1) << 1)
  307. #  define long_rsh(x,ng16) ((uint)((x) >> 16) >> (ng16 - 16))
  308. #else
  309. #  define long_lsh3(x) ((x) << 3)
  310. #  define long_rsh(x,ng16) ((x) >> ng16)
  311. #endif
  312.     return set_fmul2fixed_(pr,
  313.                    (xahi & 0xc0000000) +
  314.                     (long_lsh3(xahi) & 0x3ffffff8) +
  315.                     long_rsh(xalo, 29),
  316.                    b);
  317. }
  318.  
  319. #endif
  320.  
  321. #if USE_FPU_FIXED
  322.  
  323. /*
  324.  * Convert from floating point to fixed point with scaling.
  325.  * These are efficient algorithms for FPU-less machines.
  326.  */
  327. #define mbits_float 23
  328. #define mbits_double 20
  329. int
  330. set_float2fixed_(fixed *pr, long /*float*/ vf, int frac_bits)
  331. {    fixed mantissa;
  332.     int shift;
  333.  
  334.     if ( !(vf & 0x7f800000) )
  335.       { *pr = fixed_0;
  336.         return 0;
  337.       }
  338.     mantissa = (fixed)((vf & 0x7fffff) | 0x800000);
  339.     shift = ((vf >> 23) & 255) - (127 + 23) + frac_bits;
  340.     if ( shift >= 0 )
  341.     {    if ( shift >= sizeof(fixed) * 8 - 24 )
  342.           return_error(gs_error_limitcheck);
  343.         if ( vf < 0 )
  344.           mantissa = -mantissa;
  345.         *pr = (fixed)(mantissa << shift);
  346.     }
  347.     else
  348.       *pr = (shift < -24 ? fixed_0 :
  349.          vf < 0 ? -(fixed)(mantissa >> -shift) :  /* truncate */
  350.          (fixed)(mantissa >> -shift));
  351.     return 0;
  352. }
  353. int
  354. set_double2fixed_(fixed *pr, ulong /*double lo*/ lo,
  355.   long /*double hi*/ hi, int frac_bits)
  356. {    fixed mantissa;
  357.     int shift;
  358.  
  359.     if ( !(hi & 0x7ff00000) )
  360.       { *pr = fixed_0;
  361.         return 0;
  362.       }
  363.     /* We only use 31 bits of mantissa even if sizeof(long) > 4. */
  364.     mantissa = (fixed)(((hi & 0xfffff) << 10) | (lo >> 22) | 0x40000000);
  365.     shift = ((hi >> 20) & 2047) - (1023 + 30) + frac_bits;
  366.     if ( shift > 0 )
  367.       return_error(gs_error_limitcheck);
  368.     *pr = (shift < -30 ? fixed_0 :
  369.            hi < 0 ? -(fixed)(mantissa >> -shift) :  /* truncate */
  370.            (fixed)(mantissa >> -shift));
  371.     return 0;
  372. }
  373. /*
  374.  * Given a fixed value x with fbits bits of fraction, set v to the mantissa
  375.  * (left-justified in 32 bits) and f to the exponent word of the
  376.  * corresponding floating-point value with mbits bits of mantissa in the
  377.  * first word.  (The exponent part of f is biased by -1, because we add the
  378.  * top 1-bit of the mantissa to it.)
  379.  */
  380. static const byte f2f_shifts[] =
  381.  { 4, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0 };
  382. #define f2f_declare(v, f)\
  383.     ulong v;\
  384.     long f
  385. #define f2f(x, v, f, mbits, fbits)\
  386.     if ( x < 0 )\
  387.       f = 0xc0000000 + (29 << mbits) - ((long)fbits << mbits), v = -x;\
  388.     else\
  389.       f = 0x40000000 + (29 << mbits) - ((long)fbits << mbits), v = x;\
  390.     if ( v < 0x8000 )\
  391.       v <<= 15, f -= 15 << mbits;\
  392.     if ( v < 0x800000 )\
  393.       v <<= 8, f -= 8 << mbits;\
  394.     if ( v < 0x8000000 )\
  395.       v <<= 4, f -= 4 << mbits;\
  396.     { int shift = f2f_shifts[v >> 28];\
  397.       v <<= shift, f -= shift << mbits;\
  398.     }
  399. long
  400. fixed2float_(fixed x, int frac_bits)
  401. {    f2f_declare(v, f);
  402.  
  403.     if ( x == 0 )
  404.       return 0;
  405.     f2f(x, v, f, mbits_float, frac_bits);
  406.     return f + (((v >> 7) + 1) >> 1);
  407. }
  408. void
  409. set_fixed2double_(double *pd, fixed x, int frac_bits)
  410. {    f2f_declare(v, f);
  411.  
  412.     if ( x == 0 )
  413.       { ((long *)pd)[1 - arch_is_big_endian] = 0;
  414.         ((ulong *)pd)[arch_is_big_endian] = 0;
  415.       }
  416.     else
  417.       { f2f(x, v, f, mbits_double, frac_bits);
  418.         ((long *)pd)[1 - arch_is_big_endian] = f + (v >> 11);
  419.         ((ulong *)pd)[arch_is_big_endian] = v << 21;
  420.       }
  421. }
  422.  
  423. /*
  424.  * Compute A * B / C when 0 <= B < C and A * B exceeds (or might exceed)
  425.  * the capacity of a long.
  426.  */
  427. #ifdef DEBUG
  428. struct { long mnanb, mnab, manb, mab, mnc, mdq, mde, mds, mqh, mql; } fmq_stat;
  429. #  define mincr(x) ++fmq_stat.x
  430. #else
  431. #  define mincr(x) DO_NOTHING
  432. #endif
  433. fixed
  434. fixed_mult_quo(fixed signed_A, fixed B, fixed C)
  435. {    /* First compute A * B in double-fixed precision. */
  436.     ulong A = (signed_A < 0 ? -signed_A : signed_A);
  437.     long msw;
  438.     ulong lsw;
  439.     ulong p1;
  440.  
  441. #define num_bits (sizeof(fixed) * 8)
  442. #define half_bits (num_bits / 2)
  443. #define half_mask ((1L << half_bits) - 1)
  444.     if ( B <= half_mask )
  445.       { if ( A <= half_mask )
  446.           { fixed Q = (ulong)(A * B) / (ulong)C;
  447.  
  448.             mincr(mnanb);
  449.         return (signed_A < 0 ? -Q : Q);
  450.           }
  451.         /*
  452.          * We might still have C <= half_mask, which we can
  453.          * handle with a simpler algorithm.
  454.          */
  455.         lsw = (A & half_mask) * B;
  456.         p1 = (A >> half_bits) * B;
  457.         if ( C <= half_mask )
  458.           { ulong q0 = (p1 += lsw >> half_bits) / C;
  459.             ulong rem = ((p1 - C * q0) << half_bits) + (lsw & half_mask);
  460.         ulong Q = (q0 << half_bits) + rem / C;
  461.  
  462.         mincr(mnc);
  463.         return (signed_A < 0 ? -Q : Q);
  464.           }
  465.         msw = p1 >> half_bits;
  466.         mincr(manb);
  467.       }
  468.     else if ( A <= half_mask )
  469.       { p1 = A * (B >> half_bits);
  470.         msw = p1 >> half_bits;
  471.         lsw = A * (B & half_mask);
  472.         mincr(mnab);
  473.       }
  474.     else
  475.       { /* We have to compute all 4 products.  :-( */
  476.         ulong lo_A = A & half_mask;
  477.         ulong hi_A = A >> half_bits;
  478.         ulong lo_B = B & half_mask;
  479.         ulong hi_B = B >> half_bits;
  480.         ulong p1x = hi_A * lo_B;
  481.  
  482.         msw = hi_A * hi_B;
  483.         lsw = lo_A * lo_B;
  484.         p1 = lo_A * hi_B;
  485.         if ( p1 > max_ulong - p1x )
  486.           msw += 1L << half_bits;
  487.         p1 += p1x;
  488.         msw += p1 >> half_bits;
  489.         mincr(mab);
  490.       }
  491.     /* Finish up by adding the low half of p1 to the high half of lsw. */
  492. #if max_fixed < max_long
  493.     p1 &= half_mask;
  494. #endif
  495.     p1 <<= half_bits;
  496.     if ( p1 > max_ulong - lsw )
  497.       msw++;
  498.     lsw += p1;            
  499.     /*
  500.      * Now divide the double-length product by C.  Note that we know msw
  501.      * < C (otherwise the quotient would overflow).  Start by shifting
  502.      * (msw,lsw) and C left until C >= 1 << (num_bits - 1).
  503.      */
  504.     { ulong denom = C;
  505.       int shift = 0;
  506.  
  507. #define bits_4th (num_bits / 4)
  508.       if ( denom < 1L << (num_bits - bits_4th) )
  509.         { mincr(mdq);
  510.           denom <<= bits_4th, shift += bits_4th;
  511.         }
  512. #undef bits_4th
  513. #define bits_8th (num_bits / 8)
  514.       if ( denom < 1L << (num_bits - bits_8th) )
  515.         { mincr(mde);
  516.           denom <<= bits_8th, shift += bits_8th;
  517.         }
  518. #undef bits_8th
  519.       while ( !(denom & (1L << (num_bits - 1))) )
  520.         { mincr(mds);
  521.           denom <<= 1, ++shift;
  522.         }
  523.       msw = (msw << shift) + (lsw >> (num_bits - shift));
  524.       lsw <<= shift;
  525. #if max_fixed < max_long
  526.       lsw &= (1L << (sizeof(fixed) * 8)) - 1;
  527. #endif
  528.       /* Compute a trial upper-half quotient. */
  529.       { ulong hi_D = denom >> half_bits;
  530.         ulong lo_D = denom & half_mask;
  531.         ulong hi_Q = (ulong)msw / hi_D;
  532.         /* hi_Q might be too high by 1 or 2, but it isn't too low. */
  533.         ulong p0 = hi_Q * hi_D;
  534.         ulong p1 = hi_Q * lo_D;
  535.         ulong hi_P;
  536.  
  537.         while ( (hi_P = p0 + (p1 >> half_bits)) > msw ||
  538.             (hi_P == msw && ((p1 & half_mask) << half_bits) > lsw)
  539.           )
  540.           { /* hi_Q was too high by 1. */
  541.         --hi_Q;
  542.         p0 -= hi_D;
  543.         p1 -= lo_D;
  544.         mincr(mqh);
  545.           }
  546.         p1 = (p1 & half_mask) << half_bits;
  547.         if ( p1 > lsw )
  548.           msw--;
  549.         lsw -= p1;
  550.         msw -= hi_P;
  551.         /* Now repeat to get the lower-half quotient. */
  552.         msw = (msw << half_bits) + (lsw >> half_bits);
  553. #if max_fixed < max_long
  554.         lsw &= half_mask;
  555. #endif
  556.         lsw <<= half_bits;
  557.         { ulong lo_Q = (ulong)msw / hi_D;
  558.           long Q;
  559.  
  560.           p1 = lo_Q * lo_D;
  561.           p0 = lo_Q * hi_D;
  562.           while ( (hi_P = p0 + (p1 >> half_bits)) > msw ||
  563.               (hi_P == msw && ((p1 & half_mask) << half_bits) > lsw)
  564.             )
  565.         { /* lo_Q was too high by 1. */
  566.           --lo_Q;
  567.           p0 -= hi_D;
  568.           p1 -= lo_D;
  569.           mincr(mql);
  570.         }
  571.           Q = (hi_Q << half_bits) + lo_Q;
  572.           return (signed_A < 0 ? -Q : Q);
  573.         }
  574.       }
  575.     }
  576. #undef half_bits
  577. #undef half_mask
  578. }
  579.  
  580. #endif
  581.  
  582. /* Trace calls on sqrt when debugging. */
  583. #undef sqrt
  584. #if !defined(__amigaos__)
  585. extern double sqrt(P1(double));
  586. #endif
  587. double
  588. gs_sqrt(double x, const char *file, int line)
  589. {    if ( gs_debug_c('~') )
  590.       {    fprintf(stdout, "[~]sqrt(%g) at %s:%d\n",
  591.             x, (const char *)file, line);
  592.         fflush(stdout);
  593.       }
  594.     return sqrt(x);
  595. }
  596.  
  597. /*
  598.  * Define sine and cosine functions that take angles in degrees rather than
  599.  * radians, and that are implemented efficiently on machines with slow
  600.  * (or no) floating point.
  601.  */
  602. #if USE_FPU < 0            /****** maybe should be <= 0 ? ******/
  603.  
  604. #define sin0 0.00000000000000000
  605. #define sin1 0.01745240643728351
  606. #define sin2 0.03489949670250097
  607. #define sin3 0.05233595624294383
  608. #define sin4 0.06975647374412530
  609. #define sin5 0.08715574274765817
  610. #define sin6 0.10452846326765346
  611. #define sin7 0.12186934340514748
  612. #define sin8 0.13917310096006544
  613. #define sin9 0.15643446504023087
  614. #define sin10 0.17364817766693033
  615. #define sin11 0.19080899537654480
  616. #define sin12 0.20791169081775931
  617. #define sin13 0.22495105434386498
  618. #define sin14 0.24192189559966773
  619. #define sin15 0.25881904510252074
  620. #define sin16 0.27563735581699916
  621. #define sin17 0.29237170472273671
  622. #define sin18 0.30901699437494740
  623. #define sin19 0.32556815445715670
  624. #define sin20 0.34202014332566871
  625. #define sin21 0.35836794954530027
  626. #define sin22 0.37460659341591201
  627. #define sin23 0.39073112848927377
  628. #define sin24 0.40673664307580015
  629. #define sin25 0.42261826174069944
  630. #define sin26 0.43837114678907740
  631. #define sin27 0.45399049973954675
  632. #define sin28 0.46947156278589081
  633. #define sin29 0.48480962024633706
  634. #define sin30 0.50000000000000000
  635. #define sin31 0.51503807491005416
  636. #define sin32 0.52991926423320490
  637. #define sin33 0.54463903501502708
  638. #define sin34 0.55919290347074679
  639. #define sin35 0.57357643635104605
  640. #define sin36 0.58778525229247314
  641. #define sin37 0.60181502315204827
  642. #define sin38 0.61566147532565829
  643. #define sin39 0.62932039104983739
  644. #define sin40 0.64278760968653925
  645. #define sin41 0.65605902899050728
  646. #define sin42 0.66913060635885824
  647. #define sin43 0.68199836006249848
  648. #define sin44 0.69465837045899725
  649. #define sin45 0.70710678118654746
  650. #define sin46 0.71933980033865108
  651. #define sin47 0.73135370161917046
  652. #define sin48 0.74314482547739413
  653. #define sin49 0.75470958022277201
  654. #define sin50 0.76604444311897801
  655. #define sin51 0.77714596145697090
  656. #define sin52 0.78801075360672190
  657. #define sin53 0.79863551004729283
  658. #define sin54 0.80901699437494745
  659. #define sin55 0.81915204428899180
  660. #define sin56 0.82903757255504174
  661. #define sin57 0.83867056794542394
  662. #define sin58 0.84804809615642596
  663. #define sin59 0.85716730070211222
  664. #define sin60 0.86602540378443860
  665. #define sin61 0.87461970713939574
  666. #define sin62 0.88294759285892688
  667. #define sin63 0.89100652418836779
  668. #define sin64 0.89879404629916704
  669. #define sin65 0.90630778703664994
  670. #define sin66 0.91354545764260087
  671. #define sin67 0.92050485345244037
  672. #define sin68 0.92718385456678731
  673. #define sin69 0.93358042649720174
  674. #define sin70 0.93969262078590832
  675. #define sin71 0.94551857559931674
  676. #define sin72 0.95105651629515353
  677. #define sin73 0.95630475596303544
  678. #define sin74 0.96126169593831889
  679. #define sin75 0.96592582628906831
  680. #define sin76 0.97029572627599647
  681. #define sin77 0.97437006478523525
  682. #define sin78 0.97814760073380558
  683. #define sin79 0.98162718344766398
  684. #define sin80 0.98480775301220802
  685. #define sin81 0.98768834059513777
  686. #define sin82 0.99026806874157036
  687. #define sin83 0.99254615164132198
  688. #define sin84 0.99452189536827329
  689. #define sin85 0.99619469809174555
  690. #define sin86 0.99756405025982420
  691. #define sin87 0.99862953475457383
  692. #define sin88 0.99939082701909576
  693. #define sin89 0.99984769515639127
  694. #define sin90 1.00000000000000000
  695.  
  696. private const double sin_table[361] = {
  697.   sin0,
  698.   sin1, sin2, sin3, sin4, sin5, sin6, sin7, sin8, sin9, sin10,
  699.   sin11, sin12, sin13, sin14, sin15, sin16, sin17, sin18, sin19, sin20,
  700.   sin21, sin22, sin23, sin24, sin25, sin26, sin27, sin28, sin29, sin30,
  701.   sin31, sin32, sin33, sin34, sin35, sin36, sin37, sin38, sin39, sin40,
  702.   sin41, sin42, sin43, sin44, sin45, sin46, sin47, sin48, sin49, sin50,
  703.   sin51, sin52, sin53, sin54, sin55, sin56, sin57, sin58, sin59, sin60,
  704.   sin61, sin62, sin63, sin64, sin65, sin66, sin67, sin68, sin69, sin70,
  705.   sin71, sin72, sin73, sin74, sin75, sin76, sin77, sin78, sin79, sin80,
  706.   sin81, sin82, sin83, sin84, sin85, sin86, sin87, sin88, sin89, sin90,
  707.   sin89, sin88, sin87, sin86, sin85, sin84, sin83, sin82, sin81, sin80,
  708.   sin79, sin78, sin77, sin76, sin75, sin74, sin73, sin72, sin71, sin70,
  709.   sin69, sin68, sin67, sin66, sin65, sin64, sin63, sin62, sin61, sin60,
  710.   sin59, sin58, sin57, sin56, sin55, sin54, sin53, sin52, sin51, sin50,
  711.   sin49, sin48, sin47, sin46, sin45, sin44, sin43, sin42, sin41, sin40,
  712.   sin39, sin38, sin37, sin36, sin35, sin34, sin33, sin32, sin31, sin30,
  713.   sin29, sin28, sin27, sin26, sin25, sin24, sin23, sin22, sin21, sin20,
  714.   sin19, sin18, sin17, sin16, sin15, sin14, sin13, sin12, sin11, sin10,
  715.   sin9, sin8, sin7, sin6, sin5, sin4, sin3, sin2, sin1, sin0,
  716.   -sin1, -sin2, -sin3, -sin4, -sin5, -sin6, -sin7, -sin8, -sin9, -sin10,
  717. -sin11, -sin12, -sin13, -sin14, -sin15, -sin16, -sin17, -sin18, -sin19, -sin20,
  718. -sin21, -sin22, -sin23, -sin24, -sin25, -sin26, -sin27, -sin28, -sin29, -sin30,
  719. -sin31, -sin32, -sin33, -sin34, -sin35, -sin36, -sin37, -sin38, -sin39, -sin40,
  720. -sin41, -sin42, -sin43, -sin44, -sin45, -sin46, -sin47, -sin48, -sin49, -sin50,
  721. -sin51, -sin52, -sin53, -sin54, -sin55, -sin56, -sin57, -sin58, -sin59, -sin60,
  722. -sin61, -sin62, -sin63, -sin64, -sin65, -sin66, -sin67, -sin68, -sin69, -sin70,
  723. -sin71, -sin72, -sin73, -sin74, -sin75, -sin76, -sin77, -sin78, -sin79, -sin80,
  724. -sin81, -sin82, -sin83, -sin84, -sin85, -sin86, -sin87, -sin88, -sin89, -sin90,
  725. -sin89, -sin88, -sin87, -sin86, -sin85, -sin84, -sin83, -sin82, -sin81, -sin80,
  726. -sin79, -sin78, -sin77, -sin76, -sin75, -sin74, -sin73, -sin72, -sin71, -sin70,
  727. -sin69, -sin68, -sin67, -sin66, -sin65, -sin64, -sin63, -sin62, -sin61, -sin60,
  728. -sin59, -sin58, -sin57, -sin56, -sin55, -sin54, -sin53, -sin52, -sin51, -sin50,
  729. -sin49, -sin48, -sin47, -sin46, -sin45, -sin44, -sin43, -sin42, -sin41, -sin40,
  730. -sin39, -sin38, -sin37, -sin36, -sin35, -sin34, -sin33, -sin32, -sin31, -sin30,
  731. -sin29, -sin28, -sin27, -sin26, -sin25, -sin24, -sin23, -sin22, -sin21, -sin20,
  732. -sin19, -sin18, -sin17, -sin16, -sin15, -sin14, -sin13, -sin12, -sin11, -sin10,
  733.   -sin9, -sin8, -sin7, -sin6, -sin5, -sin4, -sin3, -sin2, -sin1, -sin0
  734. };
  735.  
  736. double
  737. gs_sin_degrees(double ang)
  738. {    int ipart;
  739.     if ( is_fneg(ang) )
  740.       ang = 180 - ang;
  741.     ipart = (int)ang;
  742.     if ( ipart >= 360 )
  743.       { int arem = ipart % 360;
  744.         ang -= (ipart - arem);
  745.         ipart = arem;
  746.       }
  747.     return
  748.       (ang == ipart ? sin_table[ipart] :
  749.        sin_table[ipart] + (sin_table[ipart+1] - sin_table[ipart]) *
  750.          (ang - ipart));
  751. }
  752.  
  753. double
  754. gs_cos_degrees(double ang)
  755. {    int ipart;
  756.     if ( is_fneg(ang) )
  757.       ang = 90 - ang;
  758.     else
  759.       ang += 90;
  760.     ipart = (int)ang;
  761.     if ( ipart >= 360 )
  762.       { int arem = ipart % 360;
  763.         ang -= (ipart - arem);
  764.         ipart = arem;
  765.       }
  766.     return
  767.       (ang == ipart ? sin_table[ipart] :
  768.        sin_table[ipart] + (sin_table[ipart+1] - sin_table[ipart]) *
  769.          (ang - ipart));
  770. }
  771.  
  772. void
  773. gs_sincos_degrees(double ang, gs_sincos_t *psincos)
  774. {    psincos->sin = gs_sin_degrees(ang);
  775.     psincos->cos = gs_cos_degrees(ang);
  776. }
  777.  
  778. #else                /* we have floating point */
  779.  
  780. static const int isincos[5] = { 0, 1, 0, -1, 0 };
  781.  
  782. double
  783. gs_sin_degrees(double ang)
  784. {    double quot = ang / 90;
  785.     if ( floor(quot) == quot )
  786.       { int quads = (int)fmod(quot, 4) & 3;    /* & 3 because might be < 0 */
  787.         return isincos[quads];
  788.       }
  789.     return sin(ang * (M_PI / 180));
  790. }
  791.  
  792. double
  793. gs_cos_degrees(double ang)
  794. {    double quot = ang / 90;
  795.     if ( floor(quot) == quot )
  796.       { int quads = (int)fmod(quot, 4) & 3;    /* & 3 because might be < 0 */
  797.         return isincos[quads + 1];
  798.       }
  799.     return cos(ang * (M_PI / 180));
  800. }
  801.  
  802. void
  803. gs_sincos_degrees(double ang, gs_sincos_t *psincos)
  804. {    double quot = ang / 90;
  805.     if ( floor(quot) == quot )
  806.       { int quads = (int)fmod(quot, 4) & 3;    /* & 3 because might be < 0 */
  807.         psincos->sin = isincos[quads];
  808.         psincos->cos = isincos[quads + 1];
  809.       }
  810.     else
  811.       { double arad = ang * (M_PI / 180);
  812.         psincos->sin = sin(arad);
  813.         psincos->cos = cos(arad);
  814.       }
  815. }
  816.  
  817. #endif                /* USE_FPU */
  818.